function VARbc = doVAR_BiasCorrected(VAR,nboot,clevel,tol,delta_update)
% Computes bias corrected RF VAR coefficients beta from bootstrap.
% Also computes Cholesky IRF from bias corrected VAR.

res       = detrend(VAR.res,'constant');
conv_test = 1;
bet_bs    = VAR.bet0(1:VAR.Ny*VAR.p,:);                                     % initial guess for beta: estimated beta

while conv_test > tol;                                                                   
    
    bet_bs_smp = zeros(VAR.Ny*VAR.p,VAR.Ny,nboot);
    
    % inner loop: wild bootstrap 
    for jj=1:nboot;
        
        %% (i) draw bootstrap sample
        rr    = 1-2*(rand(VAR.T,1)>0.5);                                    
        resb  = (res.*(rr*ones(1,VAR.Ny)));                                % draw residuals 
           t  = 1;
           Yb = zeros(VAR.T,VAR.Ny);
           Xb = zeros(VAR.T,VAR.Ny*VAR.p); 
           
        for i=1:VAR.N;                                                     % loop over each country, first row for each country set to actual value          
            % initial condition for country i
            Xb(t,:) = VAR.X(t,1:VAR.Ny*VAR.p) ;                                
            for j=t:(t+VAR.Ti(i)-1);                                       % loop over time within each country
                 Yb(j,:)   = Xb(j,:)*bet_bs(1:VAR.Ny*VAR.p,:) + resb(j,:); % follows de Vos, Evaraert, Ruyssen             
                 if j<t+VAR.Ti(i)-1;
                 Xb(j+1,:) = [Yb(j,:), Xb(j,1:((VAR.p-1)*VAR.Ny))];  
                 end
            end 
            t = t + VAR.Ti(i);
        end
        
        %% (ii) estimate VAR with this boostrappred data
        bet_tmp            = [Xb VAR.D(:,1:(VAR.N-1)) ones(length(Xb),1)] \ Yb;
        bet_bs_smp(:,:,jj) = bet_tmp(1:VAR.Ny*VAR.p,:); 
            
    end
    
    % mean beta over J bootstrap samples
    betbs_mean = mean(bet_bs_smp,3);                                       
    % compare mean beta with biased FE beta
    diff       = VAR.bet0(1:VAR.Ny*VAR.p,:) - betbs_mean(1:VAR.Ny*VAR.p,:) % (OLS - bootstrapped beta) should be ~0 if the data were generated from the unbiased beta    
    diffvec    = reshape(diff, [VAR.Ny*VAR.p*VAR.Ny , 1]); 
    display('Convergence test statistic (root mean squared difference)=')
    conv_test  = sqrt( mean(diffvec.^2 ))
    
    % update VAR coefficients (only necessary if conv_test > tolerance)
    if conv_test > tol;
    bet_bs = bet_bs + delta_update.*diff;
    end
    
end    

% define bias corrected (bc) VAR structure
%VARbc      = VAR;
VARbc.bet  = bet_bs;
%VARbc.res  = VAR.Ydm - [VAR.Xdm ]*VARbc.bet;                             % use demeaned variables 
%VARbc.Sigma= VARbc.res'*VARbc.res / (VARbc.T - VARbc.Ny*VARbc.p - VARbc.N);
%VARbc      = irfvar(VARbc);  

end
